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ABSTRACT 


RytovProp  is  a  new  approach  to  the  task  of  generating  a  large  number  of  random  realizations 
manifesting  some  aspect(s)  of  the  effect  of  turbulence  on  optical  propagation.  This  method  has 
been  applied  to  the  evaluation  of  Up-Link  performance — delivery  of  laser  power  from  a  simple 
ground  transmitter  to  a  satellite.  This  computational  method  allows  the  development  of  hundreds  of 
thousands  of  statistically  independent  random  realizations  of  the  laser  power  density  at  the  satellite 
in  just  one  or  two  minutes  on  an  ordinary  PC. 

The  RytovProp  method  is  based  on  use  of 

—  analytic  results  for  the  phase  structure  function,  the  log-amplitude  covariance,  and  the  phasedog- 
amplitude  cross-covariance,  that  have  been  developed  using  the  Rytov  approximation, 

—  the  assumption  that  turbulence  induced  phase  and  log-amplitude  perturbations  are  jointly  gaus- 
sian  random  variables,  and 

—  a  “little  trick”  from  statistics/matrix-theory  that  allows  a  realization  of  a  set  of  random  vari¬ 
ables  to  be  very  easily/quickly  developed  given  the  covariance  matrix  relating  all  of  the  random 
variables  with  each  other. 

The  “little  trick”  is  the  following.  Given  a  covariance  matrix,  and  using  the  fact  that  for  any 
covariance  matrix  it  is  possible  to  generate  a  real  matrix  which  when  multiplied  by  its  transpose 
will  be  equal  to  the  covariance  matrix — a  matrix  which  may  be  considered  to  be  the  square  root 
of  the  covariance  matrix,  then  it  can  be  shown  that  this  square  root  matrix  when  it  multiplies 
a  column  vector  of  normally  distributed,  statistically  independent  random  values  will  produce  a 
column  vector  of  gaussian  random  values  for  which  the  covariance  between  any  pair  of  elements  of 
this  column  vector  will  be  equal  to  the  corresponding  element  of  the  original  covariance  matrix. 

This  means  that  if  the  covariance  matrix  is  one  relating  the  turbulence  perturbed  phases  and  log- 
amplitudes  of  the  contributions  to  the  optical  field  at  the  satellite  from  each  of  an  array  of  points 
covering  the  transmitter  aperture,  then  these  elements  of  the  product  column  vector  (plus  any 
applicable  mean  values)  can  be  taken  as  suitably  chosen  random  realizations  of  these  turbulence 
perturbed  phases  and  log-amplitudes  of  the  contributions  to  the  optical  field  at  the  satellite  from 
each  of  an  array  of  points  covering  the  transmitter  aperture.  From  these  random  phase  and  log- 
amplitude  values  the  laser  power  density  at  the  satellite  can  be  calculated. 

With  suitable  modifications  this  computational  method  can  be  made  to  incorporate  the  effect  of 
rapid  tip/tilt  adjustment  of  the  transmitted  laser  beam — these  adjustments  being  based  on  tracking 
of  the  satellite.  Such  modifications  of  the  method  can  be  made  to  incorporate  the  effects  of  both 
anisoplanatism  and  of  finite  tip/tilt  tracking  servo  bandwidth. 

In  this  paper  the  relevant  theory  is  presented  along  with  sample  numerical  results  comparing  prob¬ 
ability  distributions  developed  using  the  RytovProp  method  and  using  the  (orders  of  magnitude 
slower)  split-step  wave  optics  propagation  method — the  comparisons  generally  showing  good  agree¬ 
ment. 


1.  Introduction 


I  use  the  term  “Up-Link”  to  refer  to  a  ground-to-space  optical  communications  system — a  system 
that  establishes  a  link  from  a  ground  based  laser  transmitter  to  a  receiver  on  an  orbiting  satellite. 
I  take  as  the  “Up-Link  Problem”  the  task  of  developing  statistical  results  characterizing  the  effects 
of  atmospheric  turbulence  on  the  signal  strength  of  such  a  link — in  particular  the  task  of  developing 
results  for  the  signal  strength’s  probability  distribution. 

The  Up-Link  Problem  is  a  subject  that  has  been  of  interest  for  many  years.  Most  of  the  prior  work 
has  been  devoted  to  attempts  to  develop  analytic  formulations — with  only  limited  success.  Because 
of  analytic  difficulties/limitations  when  quantitative  results  were  required  a  great  deal  of  reliance  had 
to  be  placed  on  a  Monte  Carlo  based  approach  with  each  random  realization  of  the  signal  strength 
being  developed  using  wave  optics  propagation  simulation  computations. 

In  preparation  for  such  wave  optics  propagation  computations 

—  the  propagation  path  is  divided  into  a  (presumably  large)  number,  N,  of  segments — with  the 
locations  along  the  propagation  path  of  N  associated  mid-segment  planes  being  defined  (along 
with  an  aperture-plane  located  at  the  ground  end  of  the  propagation  path  and  a  receiver-plane 
located  at  the  satellite  end  of  the  propagation  path,  for  a  total  of  N  +  2  planes), 

—  a  two-dimensional  grid  pattern  count-size  (for  example  1024-by-1024)  is  selected, 


and 


—  a  set  of  physical  length  that  are  to  be  associated  with  the  two  sides  of  each  of  the  N  +  2  grid 
patterns  is  defined. 

Also  in  preparation  for  wave  optics  propagation  simulation  computations 

—  a  scalar-field  representation  of  the  laser  beam  as  it  would  leave  the  aperture  of  the  ground 
based  transmitter  system  (as  it  would  leave — but  before  any  transmitter  adaptive  optics  or 
tilt  correction  is  applied)  is  defined. 

The  Monte  Carlo  method  relies  on  wave  optics  propagation  simulation  computations  to  generate  a 
large  number,  K,  of  random  realizations  of  the  scalar  optical  field  at  the  receiver-plane,  calculates 
the  signal  strength  that  is  to  be  associated  with  this  scalar  field  for  each  of  the  K  realizations,  and 
then  forms  an  estimate  of  the  statistical  parameters  that  characterize  the  random  variations  of  the 
signal  strength  from  these  K  different  realizations  of  the  signal  strength.  If  K  is  large  enough  and  if 
there  is  sufficient  statistical  independence  amongst  these  K  realizations  of  the  strength  of  the  signal, 
then  the  estimates  of  these  statistical  parameters  will  be  sufficiently  accurate. 

Regarding  the  statistical  independence  of  this  collection  of  K  results  for  the  signal  strength  the 
following  should  be  noted. 

—  The  wave  optics  propagation  calculations  are  carried  out  with  a  set  of  N  randomly  chosen 
turbulence  pattern  realizations,  one  for  each  of  the  N  mid-segment  planes — each  of  these 
N  patterns  being  chosen  in  a  way  that  is  statistically  compatible  with  the  expected  optical 
strength  of  turbulence  in  the  corresponding  one  of  the  N  segments  of  the  propagation  path. 

—  In  most  cases  it  is  necessary  to  run  the  simulation  over  a  sequence  of  instants  of  time — shifting 
the  random  turbulence  patterns  from  instant-to-instant  in  a  manner  that  represents  the  effect 
of  line-of-sight  slewing  (associated  with  tracking  the  moving  satellite)  and  the  effect  of  the 
ambient  winds — to  properly  initialize  the  operation  of  the  adaptive  optics  and/or  the  tilt 
tracking  systems. 
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—  While  wave  optics  propagation  calculations  may  be  conducted  in  such  a  way  that  for  each 
of  the  K  results  for  the  signal  strength  a  set  of  N  turbulence  patterns  unique  to  that  value 
of  the  signal  strength  will  have  been  chosen,  in  general  the  instant-to-instant  sequence  of 
operation  is  continued  for  a  long  sequence  of  instants  (after  initialization  is  completed)  and 
for  each  instant  (except  the  first  few  which  are  required  for  initialization)  a  signal  strength 
value  is  developed — all  of  which  signal  strength  values  go  into  making  up  the  set  of  K  signal 
strength  values  used  in  forming  the  estimated  values  of  the  statistical  parameters.  This  has 
the  negative  effect  that 

-  successive  values  of  the  signal  strength  are  strongly  correlated,  accordingly  reducing 
the  number  of  degrees-of-freedom  in  the  set  of  K  signal  strength  values 

from  which  it  follows  that 

-  the  accuracy  of  the  estimates  of  the  statistical  parameters  describing  the  random  vari¬ 
ations  of  the  signal  strength  is  reduced. 

That  a  long  sequence  of  instants,  with  associated  signal  strength  values,  computed  with  the  same 
(continuously  shifted)  set  of  N  turbulence  patterns,  is  used  to  provide  multiple  contributions  to 
the  set  of  K  signal  strength  values — that  it  is  used  despite  the  consequent  reduction  in  the  number 
of  degrees-of-freedom  in  the  set  of  K  signal  strength  values  and  degradation  of  the  accuracy  of 
the  estimated  statistical  parameters — reflects  the  general  slowness  of  the  wave  optics  propagation 
computational  method.  To  generate  K  statistically  independent  values  for  the  signal  strength  would 
require  considerably  more  time,  and  this  is  generally  not  done  because  of  the  relatively  large  amount 
of  time  it  takes  to  develop  wave  optics  propagation  results  for  each  instant. 

The  RytovProp  method,  the  presentation  of  which  is  the  subject  of  this  paper,  has  been  designed 
to  produce  a  large  set  of  statistically  independent  random  results  for  the  signal  strength,  and  to  be 
able  to  carry  this  out  very  rapidly  ( i.e .  orders  of  magnitude  more  rapidly  than  can  be  done  using 
a  wave  optics  propagation  based  approach),  thus  supporting  a  Monte  Carlo  based  approach  to  the 
Up-Link  Problem.  It  is  to  be  noted  that,  in  distinction  to  wave  optics  propagation  simulation,  the 
soundness  of  the  analytic  basis  for  the  RytovProp  method  is  limited  to  the  propagation  regime  for 
which  the  value  of  the  Rytov  number  (i.e.  the  value  of  the  log-amplitude  variance  computed  using 
the  Rytov  approximation)  is  significantly  less  than  unity. 


2.  Calculating  The  Signal  Strength  At  The  Receiver  Plane 

For  RytovProp  calculations  I  consider  a  square  patterned  lattice  of  points  on  the  aperture  plane  of 
the  ground  based  laser  transmitter,  an  array  of  points  sufficiently  dense  that  it  can  be  considered 
to  adequately  sample  the  aperture.  I  chose  the  spacing  between  adjacent  lattice  points  to  be  such 
that  there  are  P  lattice  points  within  the  aperture,  denoting  their  positions  by  rp  =  (x  ,y  )  where 
p  =  {1,  2,  3,  . . .  ,  P}.  I  shall  use  the  notation  r0  to  denote  the  center  of  the  aperture  and  will 
adjust/shift  the  placement  of  this  lattice  of  points  so  that  it  is  symmetric  about  this  center-point 
but  does  not  include  that  center-point — so  the  center-point  (at  r0)  will  be  in  the  center  of  the 
small  square  formed  by  the  four  points  of  the  lattice  that  are  closest  to  r0 .  There  are  thus  P  +  1 
points/positions  defined  on  the  aperture  of  the  laser  transmitter. 

I  use  the  notation  t0  to  denote  what  I  call  the  “current  time  on  the  ground.”  By  this  phrase  I  mean 
the  time  to  be  associated  with  a  value  of  the  signal  strength  at  the  receiver  at  some  instant,  but  I 
take  this  time  as  being  that  of  the  instant,  at  the  ground,  when  that  signal  was  transmitted.  (Thus 
when  I  speak  of  a  time,  t0 ,  this  time  is  to  be  associated  with  the  signal  at  the  receiver  on  the  satellite 
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but  it  is  to  be  understood  that  what  is  being  referred  to  is  the  time  when  that  received  signal  left 
the  ground  based  transmitter.)  There  is  of  course  some  time  required  for  light  to  propagate  from  the 
ground  to  the  satellite — so  the  time  when  the  signal  actually  is  received  at  the  satellite  is  latter  than 
t0,  but  that  latter  time  per  se  plays  no  part  in  the  analysis  that  I  am  presenting  here.  Associated 
with  the  current  time,  t0 ,  I  will  use  the  notation  Q0  to  denote  the  direction  in  which  the  laser  signal 
is  transmitted  at  time  t0 . 

I  shall  consider  a  sequence  of  M  prior  times,  tm ,  where  m  =  {1,2,  3,  . . .  ,  M}  -times  when  beacon- 
signals,  originating  from  the  satellite,  are  received  on  the  ground — received  from  the  corresponding 
sequence  of  direction,  0m.  It  is  to  be  noted  that  these  directions  vary  linearly  with  the  time.  To 
extend  this  understanding  regarding  the  direction  notation  so  that  it  will  also  apply  to  the  current 
time,  t0,  and  its  associated  direction,  0o,  it  is  to  be  noted  that  in  its  relation  to  all  the  other 
directions,  0m,  the  value  of  90  follows  this  same  linear  variation  with  time  except  for  the  presence 
of  an  additive  constant,  a  constant  that  is  equal  to  the  so  called  point-ahead  angle.  For  the  current 
time,  t0, 1  am  interested  in  transmission  of  the  laser  beam  to  the  satellite  rather  than  in  reception  of 
the  beacon  from  the  satellite,  so  the  corresponding  (transmission)  direction,  0Q ,  must  be  off-set  from 
the  direction  in  which  a  beacon  would  be  received  from  the  satellite  at  that  time  by  the  point-ahead 
angle — the  angle  that  a  transmitted  laser  beam  nominally  has  to  lead  the  apparent  direction  to  the 
satellite  if  the  laser  beam  is  to  be  incident  on  the  satellite. 

The  RytovProp  method  considers  optical  propagation  between  each  of  the  P  +  1  points,  rp ,  where 
p  =  {0,  1,  2,  3,  ...  ,  P}  on  the  plane  of  the  laser  transmitter  aperture  and  a  receiver/beacon  point 
on  the  satellite — which  receiver/beacon  point,  because  of  the  satellite’s  motion,  has  a  position  that 
depends  on  time,  being  different  for  each  of  the  M  + 1  times  in  the  set  consisting  of  the  current  time, 
t0,  and  the  M  prior  times,  tm  where  m  =  {  1,  2,  3,  . . .  ,  M}  —i.e.  the  set  trn  with  m  =  {0,  1,  2, 
3,  . . .  ,  M}.  I  use  the  notation  <f>  m  and  tp  m  to  denote  the  turbulence  induced  perturbation  of  the 
phase  and  of  the  log  of  the  amplitude  respectively  of  the  optical  field  in  propagation  between  the 
point  at  rp  on  the  plane  of  the  transmitter  aperture  and  the  location  of  the  receiver/beacon  point 
on  the  satellite  at  time  tm.  It  is  to  be  noted  that,  by  virtue  of  the  reciprocity  principle  in  optical 
propagation  between  two  points,  it  is  not  necessary  to  indicate  which  is  the  source  point;  the  same 
values  for  the  phase  and  log-amplitude  perturbations  apply  which  ever  one  of  the  two  points  is  the 
source  point. 

If  a  statistically  appropriate  randomly  selected  realization  of  the  set  of  P  x  (M  +  1)  values  for  the 
turbulence  induced  phase  perturbations,  <j)p  m ,  and  for  the  set  of  P  values  for  the  turbulence  induced 
log-amplitude  perturbations,  £p  0 ,  were  available  then  the  Strehl  ratio  associated  with  the  signal 
strength,  S ,  could  be  calculated  according  to  the  formula  that 

P  2 

Z  Ar  A2  eXp(£P.o  +  1  [0P.O  -  V,  ]  ) 

— - p - I - •  W 

Z4A2 

p= i 

where  Ap  represents  the  amplitude  of  the  laser  beam  at  position  rp  on  the  transmitter’s  aperture, 
A  represents  the  distance  between  adjacent  points  in  the  set  of  P  points  covering  the  transmitter’s 
aperture,  and  ip  represents  phase  shift  at  position  rp  on  the  transmitter’s  aperture  that  is  to  be 
associated  with  the  adaptive  optics  and/or  tilt  corrections. 

The  value  of  <p  may  be  considered  to  be  formed  as  a  weighted  sum  of  all  of  the  m  =  {1,  2,  3, 
. . .  ,  M}  prior  time  phase  values,  4>  m,  at  all  of  the  p  =  {1,  2,  3,  ...  ,  P}  positions  on  the  aperture. 
The  m-dependence  of  the  values  of  the  weighting  factors  (along  with  the  values,  t  ,  of  the  prior 
times)  defines  the  servo  bandwidth  of  the  adaptive  optics  and/or  tilt  tracking  control  systems. 
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For  the  Up-Link  Problem  I  restrict  attention  to  transmitter  systems  with  only  tilt  tracking  controls. 
For  such  a  system,  using  the  notation  k  to  denote  the  optical  wave  number,  where 

k  =  2  7T  /  A ,  (2) 

with  A  denoting  the  laser  (and  beacon)  wave  length,  the  value  of  ipp  may  be  considered  to  be  given 
by  the  equation 

<PP  =k(xp'dX+yp'dY),  (3) 


where 
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(7) 
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M 
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/  v  m  m  ' 

771—1 

(8) 

with  the  set  of  M  different  values  of  am  representing  the  above  mentioned  weighting  factors — a  set 
of  coefficients  that  may  be  thought  of  as  defining  the  bandwidth  of  the  beacon-tip/tilt  tracking  servo 
system.  The  quantities  and  $  represent  the  two  components  of  the  beacon  tilt,  measured 

at  the  mTH  prior  time,  trn ,  and  may  be  considered  to  have  values  given  by  the  equation 

0*  = 

P 

=  fc_1  E  *P 

P=  1 

,  and 

P 

E  = fc_1  E  yp  ^P,m  • 

p=i 

(9) 

At  this  point  it  is  appropriate  to  note  that  by  virtue  of  Eq.  (6)  I  can  recast  Eq.  (9)  as 

X 

P 

= fc_1  E 

p= i 

®P  {  <AP,m  - 

}  and 

p 

E  =  ^  E  -  ^0,m  } 

P=1 

p 

=  k~l  E 

p= i 

®p  Em  > 

P 

=  E  yP  Em . 

p=i 

(10) 

where 

Em  =  ^p,m  -  K 

m  ’ 

(11) 

I  will  refer  to  the  quantity  denoted  by  cf)p  m  as  the  “adjusted  phase  perturbation.” 
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It  is  also  to  be  noted  here  that  since  exp  (i<j)0m)  |2  =  1  then  Eq.  (1)  can  be  rewritten 


as 


5  = 


Jf 

E  ap  a2  exp(  [{4o  +  i[{<t>P,  o  -¥>*]) 


P= 1 


E  A 


p= i 


E  \  exp  ( [I, +I  ]  +  *  [^,o  -  vP  ] ) 


p=i 


P  1 2 

EA 

p=i 


(12) 


where 


(13) 


with  £  denoting  the  mean  ( i.e .  the  ensemble  average)  value  of  the  turbulence  induced  log-amplitude 
perturbations,  which  mean  value  is  independent  of  the  position,  r  ,  and  of  the  time,  tm ,  and  is 
defined  by  the  equation 


'  =  <<■->• 


(14) 


I  will  refer  to  the  quantity  denoted  by  £  as  the  “adjusted  log-amplitude  perturbation,”  and  will 
refer  to  the  quantity  denoted  by  £  as  the  “log-amplitude  expected  value.” 

It  can  be  seen  that  given  a  set  of  statistically  appropriate  randomly  selected1  value  for  <j>  m  for  all 
p  =  {1,  2,  3,  . . .  ,  P}  and  all  m  =  {0,  1,  2,  3,  ...  ,  M}  -a  total  of  P  x  (M  + 1)  random  values — and 
for  £  0  for  all  p  =  {1,  2,  3,  ...  ,  P}  a  total  of  an  additional  P  random  values,  then 


—  first,  making  use  of  Eq.’s  (4)  and  (9),  correspondingly  statistically  appropriate  random  pairs 
of  value  can  be  developed  for  each  of  the  prior  time  tilts,  $  and  , 


—  next,  making  use  of  Eq.  (8),  a  correspondingly  statistically  appropriate  pair  of  values  can  be 
developed  for  the  time  weighted  average  (servo  correction)  tilts,  i9  and  i9  , 


—  following  which,  making  use  of  Eq.’s  (3)  and  (4),  the  set  of  tilt-tracking  servo  correction 
induced  phase  adjustments,  <p  ,  for  each  of  the  P  positions,  r  ,  can  be  calculated, 

Using  Eq.  (12)  I  can  then  calculate  the  value  of  the  associated  (and  statistically  appropriate)  Strehl 
ratio,  S. 


This  would  seem  to  require  the  initial  generation  of  a  total  of  P  x  (M  +  1)  +  P  =  P  x  (M  +  2) 
random  values  for  each  Strehl  ratio  value  that  is  produced.  As  will  be  developed  latter  we  can 
“short  circuit”  some  of  this  process  and  directly  generate  a  set  of  random  values  for  (j>  0 ,  £p0  for 
p  =  {1,  2,  3,  . . .  ,  P},  along  with  random  values  for  and  for  m  =  {1,  2,  3,  ...  ,  M},  a 


1  By  the  term  “statistically  appropriate  randomly  selected”  I  imply  that  the  method  used  for  choosing  values  for  0 

and  for  £p  0  results  in  values  whose  statistics  are  in  accordance  with  all  that  we  know  of  the  statistics  for  these 
random  variables.  As  will  be  developed  latter,  what  we  consider  that  we  know  of  the  statistics  is  1)  that  these  random 
variables,  taken  together,  have  a  jointly  gaussian  distribution,  2)  that  they  each  have  a  mean  value  of  zero,  3)  what 
the  value  of  the  covariance  of  „  with  d>  ,  ,  is,  4)  what  the  value  of  the  covariance  of  £n  n  with  £  ,  is,  and  what 
the  cross-covariance  between  d>  and  i  .  is. 

rp,m  v/  0 
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set  for  which  the  combination  of  random  values  will  all  be  statistically  appropriate — a  set  of  only 
P  +  P  +  M  +  M  =  2  (P  +  M)  random  values — from  which  a  corresponding,  statistically  appropriate, 
Strehl  ratio  value  can  be  developed.  This  can  be  done  directly,  i.e.  without  having  to  first  generate 
the  much  larger  set  of  randomly  selected  values  of  <j>  m  for  all  combinations  of  p  =  {1,  2,  3,  P} 
and  to  =  {1,  2,  3,  ,  M}.  But  this  will  be  covered  in  Section  4.  For  now  I  concern  myself  with 

setting  forth  the  method  of  generating  the  full  set  of  value  of  <j>p  m  along  with  the  set  of  values  of 
l  0  for  all  values  of  p  and  m.  I  take  this  up  in  the  next  section. 


3.  Generating  A  Full  Set  Of  Randomly  Selected  Statistically  Appropriate  Phase  And 
Log- Amplitude  Perturbation  Values 

I  separate  this  section  into  three  sub  sections.  In  the  first  sub  section  I  present  the  basic  covariance 
matrix  based  approach  to  the  generation  of  a  set  of  statistically  appropriate  randomly  selected  values 
for  a  set  of  zero  mean,  gaussian  random  variables — providing  that  each  of  these  random  variables 
has  a  finite  variance. 

In  the  second  sub  section  I  present  the  Rytov  approximation  based  second  moment  results  for 
turbulence  induced  phase  and  log-amplitude  perturbations,  (j>  m  and  l  m  — the  second  moment 
results  referred  to  as  the  log-amplitude  covariance,  the  phase  structure  function,  and  the  phase  :  log- 
amplitude  cross-covariance  function. 

In  the  third  sub  section  I  show  how  the  second  moment  results  for  the  turbulence  induced  phase 
and  log-amplitude  perturbations,  <j>  m  and  t  m ,  that  are  presented  in  the  second  sub  section  may  be 
used  to  develop  the  covariance  and_cross-covariance  values  for  the  adjusted  phase  and  adjusted  log- 
amplitude  perturbations,  <fip  m  and  t  m .  (It  is  to  be  noted  that  these  covariance  and  cross-covariance 
values  constitute  the  elements  of  a  covariance  matrix  that  could  be  used — in  the  way  set  forth  in 
the  first  sub  section — in  the  generation  of  a  statistically  appropriate  realization  of _a  set  of  random 
values  for  the  adjusted  phase  and  adjusted  log-amplitude  perturbations,  <j)p  m  and  ip  m.) 


3.1.  Generating  A  Set  Of  Statistically  Appropriate  Zero-Mean,  Jointly-Gaussian,  Ran¬ 
dom  Variables 

Let  zn  denote  a  set  of  n  =  {1,  2,  3,  ...  ,  N}  zero-mean,  jointly-gaussian,  random  variables,  and 
let  z  denote  a  column  vector  whose  ?rTH  element  is  zn.  Using  the  angle-bracket  notation,  (...),  to 
indicate  an  ensemble  average  the  covariance  matrix,  Cz,  for  this  set  of  random  variables  can  be  seen 
to  be  given  by  the  equation 

=  (  z  zT  )  .  (15) 

This  is  an  N-by-N  size  matrix. 

Let  U_  denote  a  matrix  whose  columns  correspond  to  the  eigen-vectors  of  Cz,  and  let  S_  denote  a 
diagonal  matrix  whose  diagonal  elements  are  the  corresponding  eigen- values  of  C_ .  Both  Uz  and 
Sz  are,  like  Cz,  N-by-N  size  matrices.  I  can  write 

cz  uz  =  Uz  Sz .  (16) 

Because  the  eigen  vectors  are  ortho-normal  Uz  is  unitary/orthogonal  and  I  can  write 

Uj  Uz  =  I,  as  well  as  Uz  uj  =  I ,  (17) 
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where  I  is  the  identity  matrix  of  size  IV-by-iV  — a  diagonal  matrix  all  of  whose  diagonal  elements 
are  equal  to  unity. 

Let  S_  denote  a  diagonal  matrix  whose  diagonal  elements  are  each  equal  to  the  square  root  of  the 
corresponding  one  of  the  diagonal  elements  of  S_,  so  that 

s,  §;  =  ss .  (is) 

Now  form  the  matrix  Cz  according  to  the  equation 

C,  =  U  A  .  (19) 

It  is  easy  to  shown  that  this  matrix,  Cz,  is  the  square  root  of  the  matrix  C_  — the  square  root  of  Cz 
in  the  sense  that  Cz  CL  =  Cz .  To  prove  this — successively  making  use  of  Eq.  (19),  (18),  and  (17) 
— I  write 

C.  Cz  =  (u,  sj  (u,  s  J  =  u.  s,  sj  uj  =  u  sz  uz  =  c„  u  uj  =  Cz .  (20) 


Let  7  be  a  column  vector  of  length  N  of  zero-mean,  unity-variance,  statistically-independent,  gaus- 
sian  random  variables,  so  that 


( 7 )  =  0 ,  and  ( 7 1  )  =  I ,  (21) 

where  0  is  a  column  vector  of  length  N  all  of  whose  elements  are  equal  to  zero,  and  as  above  I  is 
the  identity  matrix  of  size  N-by-N. 

I  assert  that  the  column  vector  £  formed  according  to  the  equation 

C  =  Czl,  (22) 

is  a  statistically  appropriate  realization  of  the  random  variable  z.  To  prove  that  this  is  so  I  first  note 
that  since  the  elements  of  7  taken  as  a  set  of  random  variables  constitute  a  realizations  of  a  jointly 
gaussian  set  or  random  variables,  and  since  the  elements  of  £  are  formed  as  variously  weighted  sums 
of  the  elements  of  7,  then  the  elements  of  C  constitute  a  realization  of  a  set  of  jointly  gaussian 
random  variables.  I  next  note  that  since  (7)  =0,  then  since  C_  is  non  random  then 

<C>  =  (C  7>  =  C,<7>  =  CI0  =  0.  (23) 


Finally,  I  note  that 

<  CC  )  =  (  (C  7)  (c,  7)T  }  =  (  Cz  77T  cj  }  =  cz  <77T  >  cj  =  CzICz  =  Cz  Cz  =  Cz  .  (24) 

These  three  just  noted  fact  indicate  that  the  values  of  the  elements  of  C,  calculated  in  accordance  with 
Eq.  (22),  are  a  statistically  appropriate  realization  of  the  random  variable  elements  of  the  column 
vector  z.  They  have  the  correct  mean  value,  the  correct  set  of  covariances  and  cross-covariances, 
and  are  jointly  gaussian. 


3.2.  Rytov  Approximation  Based  Statistics  For  Turbulence  Induced  Phase  And  Log- 
Amplitude  Perturbations 
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To  develop  values  for  the  elements  of  the  covariance  matrix  C_  when  the  elements  of  z  correspond 
to  the  various  adjusted  phase  and  adjusted  log-amplitude  values,  (f>pm  and  £p  m,  that  will  have  to  be 
available  to  allow  generation  of  Strehl  ratio  values  utilizing  to  Eq.  (12),  I  start  by  defining  the  log- 
amplitude  covariance  function,  Cu,  the  phase  structure  function,  ,  and  the  phasedog-amplitude 
cross-covariance,  C^,  for  the  various  combinations  of  aperture  plane  position,  rp  and  r  , ,  and  various 
combinations  of  times,  tm  and  t  ,  and  associated  propagation  directions,  9m  and  6  , .  These  three 
statistical  second  moment  quantities  are  defined  by  the  equations 

Cu ip,  m\ p' ,  m!)  =  (  [iPtm  -  I]  [£p,tm  -  I]  (25a) 

ip,  m;  p,  m!)  =  (  [  0p,m  -  <£p,jm,  ] 2  )  ,  (25b) 

C*tip,m\p',rri)  =  ((j)pm  -£]  Y  (25c) 

It  is  to  be  noted  that  these  quantities,  Cu(p,m',p',m'),  (p,  m;  p',  m'),  and  C4te(p,m;p,,m'),  are 

defined  in  terms  of  the  turbulence  induced  phase  and  log-amplitude  perturbations,  <f>  m  and  £  m , 

and  not  in  terms  of  the  adjusted  phase  and  adjusted  log-amplitude  perturbations,  <j)  m  and  ip  m. 

Making  use  of  the  Rytov  approximation  in  developing  a  solution  for  the  wave  propagation  equation, 
assuming  that  the  atmospheric  turbulence  pattern  is  in  accordance  with  the  Kolmogorov  theory  of 
turbulence  in  the  inertial  sub  range,  and  farther  assuming  that  the  time  dependence  of  the  turbulence 
pattern  is  governed  by  the  Taylor  theory  of  frozen  turbulence,  the  following  results  can  be  developed 
for  Cu ,  ,  and  C ^ ,  namely  that 


.  .  8.16  2 
Cu(p,m;p  ,m  )  =  k 

[  dzC^iz )  [z(l  -  z/Z)/k}5/6  F(Q(p,m;p',m';z)) , 

Jo 

(26a) 

,  .  8.16  2 
iP,  m;p,m)=  ^  ^  k 

f  dzCliz)  [z(l-  z/Z)/k]5/6  G(Q(p,m;p',m';z))  , 

Jo 

(26b) 

C^e  ip,  m\  p',  m!)  =  4?]_  k2 

[  dzCl{z)  [z(  1  -  z/Z)/k\ 5/6  H (Q(p,  to;  p',  m';  z))  , 
Jo 

(26c) 

where  the  notation  Z  denotes  the  range  from  the  ground  to  the  satellite  (i.e.  from  the  laser  transmit¬ 
ter’s  aperture  plane  to  the  laser  receiver  or  beacon  source  on  the  satellite),  the  variable  of  integration, 
z,  can  be  considered  to  denote  position  along  the  propagation  path  (from  z  =  0  at  the  ground  to 
z  =  Z  at  the  satellite),  and  C^(z)  denotes  the  value  of  the  refractive  index  structure  constant  (which 
is  a  measure  of  the  optical  strength  of  turbulence)  at  the  position  z  along  the  propagation  path. 
The  three  functions  F(Q),  G(<5),  and  H(Q)  are  defined  by  the  equations 


HQ)  = 

POO 

/  d,K  J0  (k  Q)  1— cos(k2)  , 

Jo 

(27a) 

G{Q)  = 

/  dKK~8/3  1  — J  0(kQ)  1  +  cos(av2)  , 

Jo  -1  L 

(27b) 

H{Q)  = 

POO 

/  dn  k~8^3  J0  (k  Q )  sin(Av2) . 

Jo 

(27c) 

The  quantity  Q(p,  m;  p',  m! \  z)  appearing  in  Eq.  (26)  has  a  value  given  by  the  equation 


where  V(z)  denotes  the  projection  onto  a  plane  parallel  to  the  aperture  plane  of  the  vector  repre¬ 
senting  the  wind  velocity  that  exists  at  the  range  z,  where  0m  and  6  ,  denote  directions  associated 
with  arrival  of  the  beacon’s  light  at  time  tm  and  t  ,  when  m  and/or  m!  =  {1,  2,  3,  . . .  M }  and  the 
direction  associated  with  transmission  to  the  satellite  at  time  t0  when  m  and/or  m!  =  0,  and  where 
the  quantity  £p ,  which  may  be  considered  to  be  a  sort  of  Fresnel  length,  has  a  value  given  by  the 
equation 


= 


kZ 


z  (Z  —  z) 


(29) 


An  interesting  physical  interpretation  can  be  placed  on  the  definition  of  Q{p,m\p'  ,m!  \  z)  given  by 
Eq.  (28)  — an  interpretation  as  a  separation  length  divided  by  the  normalization  length  CF.  To 
explain  what  I  am  referring  to  when  I  speak  of  a  “separation  length”  I  call  attention  to  a  pair  of  rays 
to/from  the  aperture  plane.  Ray-1  is  associated  with  (p,m),  i.e.  with  the  aperture  plane  position 
r  ,  the  time  tm,  and  the  direction  0  m.  Ray-2  is  associated  with  ( p',m '),  i.e.  with  the  aperture 
plane  position  r  , ,  the  time  tm, ,  and  the  direction  6m,  .  In  accordance  with  the  Taylor  hypothesis  of 
frozen  flow,  we  may  consider  the  turbulence  pattern  per  se  to  be  unchanging  with  time  but  to  be 
transported  by  the  ambient  wind.  At  time  tm  Ray-1  pierces  the  turbulence  pattern  that  exists  at 
the  range  z  —pierces  it  at  some  position.  At  time  t  ,  Ray-2  pierces  this  same  turbulence  pattern 
at  some  other  position.  It  is  then  appropriate  to  ask  what  is  the  separation  between  where  the  first 
ray  pierced  the  turbulence  pattern  and  where  the  second  ray  pierced  that  pattern. 


It  can  be  seen  that  the  rays  pierce  a  fixed  plane  at  range  z  at  positions  that  are  separated  by  (r  — 
r  )  (1  —z/Z)  +  (0m  —  0m, )  z.  But  since  the  turbulence  pattern  has  moved  a  distance  V(z)  (tm  —  tm, ) 
in  the  time  interval  between  the  two  piercings,  the  separation  between  where  the  turbulence  pattern 
is  pierced  by  these  two  rays  is  reduced  by  V(z)  (tm  —  t  ,).  Accordingly  the  separation  of  the  two 
piercings  of  the  turbulence  pattern  per  se  is  given  by  the  expression  (rp  —  r  ,)  (1  —  z/Z)  +  (8rn  — 
0  , )  2  —  V(z)  (tm  —  t  , )  — which  is  the  numerator  of  the  fraction  on  the  right-hand-side  of  Eq.  (28) . 


Expressing  the  trigonometric  functions  in  Eq.  (27)  in  terms  of  a  sum/difference  of  exponential 
functions  and  making  use  making  use  of  well  known  results  for  certain  definite  integrals2  it  can  be 
shown  that  the  functions  F(Q ),  G(Q),  and  H(Q ),  which  are  defined  by  Eq.  (26),  have  values  given 
by  the  equations 


F(Q) 

G(Q) 

H(Q) 


i  r(  - 1)  /I 


(i  Q2)5/6  -  3  r(  -  f )  »{  exp  (  &  * 0  lFl  (  -  1;  i  Q2  i)  } 


1  r(  e)  (\  ^,2\5/6 


( 1  Qif*  -  I  r(  -  |)  exp  (A  , i)  lFl (  _  1;  1  Q2  i)  }  _  cos  ( &  tt) 


2  r(¥) 

ir(-|)3?{  exp(ii7T?;)  iFi(-  |;l;iQ2?;)  |. 


(30a) 

i  (30b) 
(30c) 


The  hyper  geometric  functions  appearing  in  this  result  can  be  evaluated  using  the  standard  power 
series  formulation — so  long  as  Q  is  not  too  large.  With  16-digit  computational  accuracy  quite 
accurate  results  can  be  developed  for  values  of  Q  as  large  as  Q  =  10,  but  for  values  of  Q  much  larger 
than  about  Q  =  12  the  results  obtained  with  16-digit  computational  accuracy  are  very  clearly  in 
error.  To  circumvent  this  limitation/difficulty  I  have  developed  the  asymptotic  series  results  that 


F(Q)  = 


r(j) 

r(-s) 


(lQT7/6-HlQ2) 


2\  —11/6 


1  - 


8  n 

3  3 


iQ2) 


-2 


8  11  14  17  (1  r)2\~4 
3  3  3  3  V4  **  ) 


sin  (i  Q 2) 


2  I.S.  Gradshteyn  and  I.M.  Ryzhik,  Tables  of  Integrals,  Series,  and  Products ,  4th  edition  (Academic  1965  New  York); 
Eq.  (6.561.14),  p.  684  and  Eq.  (6.631.1),  p.  716 
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cos  Q2) 


(31a) 


G(Q) 

H(Q) 


3  (3  <3 


2N-17/6 


8  11  14  n  n2\~2 

3  3  3  V4  ^  ) 


8  11  14  17  20  TO,  -q2\“4 
3  3  3  3  3  V  4  / 


r(~l) 

r(f) 


(\Q2f'& +F{Q)  +  \T{ 


I)  C0S(lH 


(31b) 


1  £(i) 
5  r(i) 


-1/6 


U\Q2) 


2\  —11/6 


1  - 


8  11 
3  3 


(1Q2) 


-2 


8  11  14  17  (1  r)2\~4 
3  3  3  3  V4  ^  J 


:iQ2) 


-17/6 


8  n  14  (i  r>2\~2  _i_  8  il  M  II 20  II  n2l~4 

3  3  3  V4  ^  I  3  3  3  3  3  V4^  ) 


sin(iQ2)  , 


cos(±Q2) 


(31c) 


which  results3  very  smoothly  joint  the  results  given  by  Eq.  (30)  at  Q  =  10 . 

Before  closing  this  sub  section  I  introduce  one  additional  fact.  Based  on  the  presumption  that 
the  turbulence  induced  log-amplitude  perturbations  are  normally  distributed  it  can  be  shown  by 
consideration  of  conservation  of  energy  requirements,  that 


l=-Cu  ( p ,  to; p,  to )  =  -U)  ,  (32) 

where  the  notation  Vp  is  used  to  denote  what  has  come  to  be  called  the  Rytov-number.  4>  5 

Having  now  provided  all  the  computational  formulations  needed  to  develop  Rytov  approximation 
based  second  moment  propagation  statistics  for  the  turbulence  induced  random  phase  and  log- 
amplitude  perturbations,  (j)p  m  and  £p  0 , 1  turn  in  the  last  of  these  three  sub  sections  to  the  matter  of 
applying  these  formulations  in  evaluation  of  the  covariances  and  cross-covariances  for  the  adjusted 
phase  and  adjusted  log-amplitude  perturbations,  </>p  m  and  ip  0. 


3.3.  Using  The  Rytov  Approximation  Results  In  Generating  Covariance  And  Cross- 
Covariance  Results  For  The  Adjusted  Phase  And  Adjusted  Log-Amplitude  Perturba¬ 
tions 

It  is  obvious  that  this  quantity,  </>  ,  has  a  zero  mean  value.  The  covariance  for  this  quantity, 

C_-(p,  m\p',m!),  may  accordingly  be  considered  to  be  defined  by  the  equation 

Cu  ip,  m;  p\  to')  =  <  4>p,  m,  )  .  (33) 

An  expression  allowing  evaluation  of  the  value  of  this  quantity  can  be  easily  developed  using  the 
simple  algebraic  relationship  that  (a  —  b)  (c  —  d)  =  \  [  —  (a  —  c)2  +  (a  —  d)2  +  (b  —  c)2  —  {b  —  d)2  ] . 
It  follows  from  this  and  the  definition  of  the  phase  structure  function,  V ^ ,  given  by  Eq.  (25b)  that 

iP,  "i;  P',  to')  =  \  [  -  ip ,  to;  p',  m!)  +  ip ,  to;  0,  m!) 


+  (0:  to;  p',  to')  -  (0,  to;  0,  to')  . 


(34) 


3  The  details  of  the  derivation  of  these  asymptotic  series  results  would  take  up  too  much  space  to  be  worth  presenting 
here.  Details  will  be  provided  on  request  as  my  Tech-Note  TN-205. 

4  D.L.  Fried,  “Scaling  Laws  for  Propagation  through  Turbulence,”  Atmos.  Oceanic  Opt.  11  982-990  (1998) 

5  I  use  the  notation  73/  here,  in  place  of  the  more  customary  notation  (or  cr^),  called  the  log-amplitude  variance,  to 
take  account  of  the  fact  that  under  propagation  conditions  for  which  7 Z  is  calculated  to  have  a  value  larger  than  about 
7 Z  =  0.25  Np2  there  occurs  what  is  called  “saturation  of  scintillation”  and  the  actual  variance  of  the  log-amplitude 
does  not  increase  with  increasing  strength  of  turbulence. 
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Since  d>  „  and  i  ,  ,  each  have  mean  values  of  zero  their  cross-covariance,  C --(»,  m;»' ,  to')  can  be 
considered  to  be  defined  by  the  equation 

C}l{p,  to;  p',  to')  =  (  4>Ptm  Ip,m,  )  .  (35) 

Taking  note  of  Eq.’s  (11),  (13),  and  (25c)  it  can  be  seen  that 

cfl(p,  to;  P\  m ')  =  ip ,  to;  p',  to')  -  C*  (0,  to;  p',  to')  .  (36) 


Since  the  mean  value  of  £p  m  is  zero  the  covariance  for  a  pair  of  adjusted  log-amplitude  perturbation 
values,  for  say  tp  m  and  t  ,  m,  ,  which  covariance  I  shall  denote  by  the  notation  CH(p,  m;  p',  to'),  may 
be  considered  to  be  defined  by  the  equation 

Cu(p,  to;  p',  to')  =  ( m  7  ,  m,  )  .  (37) 

Taking  note  of  Eq.’s  (13)  and  (25a)  I  see  that  I  can  write 

Cu  (p,  to;  p',  to')  =  Cu  (p,  to;  p',  to')  .  (38) 

With  Eq.’s  (34),  (36),  and  (38)  along  with  Eq.  (32)  in  hand  I  have  expressions  for  evaluation  of 
all  the  elements  of  the  covariance  matrix,  C. ,  when  the  elements  of  z  correspond  to  the  adjusted 
phase  and  adjusted  log-amplitude  perturbations,  4>pm  and  £pm.  With  such  a  covariance  matrix  I 
can  generate  statistically  appropriate  random  realizations  for  the  values  of  (j>  m  and  t  m  — and  from 
this  can  generate  statistically  appropriate  random  realizations  of  the  Up-Link  Strehl  ratio,  S. 

Unfortunately,  because  many  prior  times,  tm  for  to.  =  {l,  2,  3,. . .  ,  M}  are  generally  required  to 
allow  proper  simulation  of  the  tilt  tracking  process  that  is  represented  by  Eq.  (8),  the  size  of  the 
required  covariance  matrix,  Cz ,  can  be  awkwardly  large.  To  avoid  this  difficulty,  in  the  next  section 
I  shall  consider  the  possibility  of  making  the  beacon  tilts  at  time  tm,  namely  and  i?  ,  elements 
of  the  column  vector  z  and  developing  statistically  appropriate  randomly  selected  values  for  these 
quantities  directly  from  application  of  Eq.  (22).  This  would  allow  me  to  drop  the  task  of  generating 
P  x  M  of  the  prior  time  adjusted  phase  perturbation,  (j>  m,  random  values,  dropping  these  elements 
from  the  column  vector  z  —generating  in  their  place  only  2  M  random  tilt  values  values,  $X  and  i)Y 
— thus  greatly  reducing  the  length  of  z  and  the  size  of  its  covariance  matrix,  C_ . 


4.  Directly  Generating  Prior  Time  Beacon  Tilt  Random  Values 

If  the  i9  and  i?  random  quantities  are  to  be  elements  of  z  then  I  will  need  expressions  allowing 
evaluation  of  the  following  covariances  and  cross-covariances,  namely 

Cxx  (to;  to')  =  (  ^  ,  )  ,  (39a) 

CXY  (to;  to')  =  (  )  ,  (39b) 

CYY  (to;  to')  =  (  )  ,  (39c) 
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c,x  (p>  m; "0  =  ( 4>P,m  E ) > 

C,Y(p,m;m,)  =  (0J>ira^#>, 

c&(p,m;m')  = 

C|Y  (p,  m;  to')  =  ( 7p  m  E  >  . 


(39d) 

(39e) 

(39f) 

(39g) 


Substituting  from  Eq.  (10)  into  Eq.  (39a),  making  a  double  sum  of  the  product  of  two  sums  (using  p 
as  the  summation  index  for  one  of  the  sums  and  p'  as  the  summation  index  for  the  other  summation) , 
then  interchanging  the  order  of  summation  and  of  ensemble  averaging,  and  finally  making  use  of 
Eq.  (33)  and  then  of  Eq.  (34),  I  can  write  for  Cxx(?n;  in')  that 


Cxx(m;  m')  =  k  2  (  ^  xp  <t>pm  ^  xp,  <£p/  m,  \ 

\  P= l  p'=  l  / 


=  *-a(  E 

\  p,p'—i 


Xn  X  ,  <Z>  m  (D  ,  , 

P  p  ’P,xn  i  p’ ,rn’ 


=  k 2  x„  x  ,  ( <j>„  m  d>  ,  ,) 

/  j  P  p'  \  ~  P,m  t  pi  m'  / 

P,P'  =  1 
P 

=  k~2  E  *p>  C4>4>  (?>  ?Tl;  TO0 

P,P'  =  1 
P 

=  E  ^  I  ~  v **  (p>  m;  p'»  m')  +  (p>  m;  o, "0 

p>p'=i 

+  2^  (0:  m;  p',  rri)  -  (0,  m;  0,  m)  }  . 

Proceeding  similarly  for  CXY(m;  m')  and  CYY(m;  m')  I  obtain  the  results  that 

p 

cXY  (m;  m')  =  53  yP>  {  -  v **  (p>  TO;  p'»  m')  +  (p>  m;  o, m') 


p,p'=i 


and 


CYY  (to;  to')  =  fc  2  53  pp  2/p,  {  -  2?^  (p,  m;  p',  to')  +  V  ^  (p,  in;  0,  to') 
Pip'=i 

+  2^  (0,  m;  p',  to')  -  2?^  (0,  to;  0,  to)  }  . 


(40) 


+  2^  (0,  to;  p',  to')  -  2?^  (0,  to;  0,  m)  }  ,  (41) 


(42) 


Following  the  same  general  sort  of  analytic  development  as  was  used  in  developing  Eq.  (40)  I  can 
write  for  C  (p,  to;  in')  that 


C^x(p,  to;  to')  =  fc  1  ^  <^Pim  tp’.m'  ^ 
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=  k  1  ^  Z  ^P.m  t'.m'  ^ 

P 

=  fc_1  Z  ^p'.m'  ) 

P/  =  l 
P 

=  k~ 1  Z  Xp'C4,t(P’rn'’P'’rn') 

P'  =  1 
P 

=  k~1  Z  *„/  {  -^.(P’TO^toO  +  P^^tojO,™') 

P,  =  l 

+  ^  (0,  m;  p',  to')  -  2>^  (0,  to;  0,  to)  }  .  (43) 

The  corresponding  result  for  C-Y(p,  m;  m')  is 

p 

C4,y  (p>  TO;  m0  = fc_1  Z  yP'  i  -  ( p > m;  p,> m/)  +  (P’ TO;  °’  m0 

p'=l 

+  (0,  to;  p',  to')  -  X>^  (0,  to;  0,  to)  }  .  (44) 

Proceeding  in  essentially  the  same  way  for  the  evaluation  of  C_  (p,  to;  to'),  only  this  time  making 
use  of  Eq.  (35)  rather  than  Eq.  (33)  and  Eq.  (36)  rather  than  Eq.  (34),  I  write 

C*-X(P>  to;  to')  =  fc"1  (  7p  m  Z  Z 

\  p'=i 

=  fc_1  ^  Z  *p'  ?p,m  Z™'  ^ 

=  fc_1  Z 

p,=l 

p 

=  k -1  Z  £,/  (p', to'; p, to) 

P'  =  l 
P 

=  k~x  Z  *„/  {  c<n  (p'.  to';  p,  to)  -  (0,  to';  p,  to)  }  .  (45) 

P'  =  l 

The  corresponding  result  for  C-Y(p,  m;  m')  is 

p 

ciY  ip,  to;  to')  =  fc_1  Z  PP'  {  C<h  (p'>  to';  P,  to)  -  C*  (0,  to';  p,  m)  }  .  (46) 

p'= i 

With  these  results  in  hand  I  can  now  consider  a  version  of  z,  the  column  vector  of  random  variables, 
that  consist  of 
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—  of  P  values  for  the  current  time  phase,  <f>  0 , 

—  of  P  values  for  the  current  time  log-amplitude,  £p  0 , 

—  of  M  values  for  the  prior  times  x-components  of  tilt,  $  ,  and 

Y 

—  of  M  values  for  the  prior  times  y-components  of  tilt,  i?  , 

a  total  of  2  (P  +  M)  random  values  needed  to  generate  a  Strehl  ratio  value  when  the  prior  time  tilts 
are  generated  directly  rather  than  calculated  from  all  of  the  prior  time  adjusted  phases  values.  This 
number,  2  (P  +  M),  is  considerably  less  than  the  P  ( M  +  2)  number  of  elements —  P  (1  +  M)  for 
the  current  time  and  all  the  prior  time  adjusted  phase  perturbation  values  and  P  for  the  current 
time  adjusted  log-amplitude  perturbation  values — that  would  have  had  to  have  been  generated  if 
the  prior  time  tilts  were  calculated  from  directly  generated  prior  time  adjusted  phase  perturbation 
values.  In  terms  of  the  covariance  matrix  and  the  computational  time  charge  for  extraction  of  its 
eigen  vectors  and  eigen  values  this  is  a  very  important  reduction  in  matrix  size.  It  is  also  a  significant 
factor  in  terms  of  the  time  required  for  generation  of  each  random  realization  of  the  Strehl  ratio’s 
value,  particularly  if  millions  of  random  realizations  are  required — as  in  the  development  of  results 
relating  to  low  probability  of  occurrence  Strehl  ratio  values. 

With  all  of  the  needed  computational  tools  defining  RytovProp  now  in  place  I  turn  in  the  next 
section  to  a  presentation  of  results  demonstrating  the  soundness  of  the  RytovProp  method.  I  will 
present  sample  results  generated  using  this  method  and  compare  these  results  with  results  generated 
utilizing  the  wave  optics  propagation  simulation  methods. 


5.  Testing/Validation  Of  The  RytovProp  Method 

To  allow  an  assessment  of  the  soundness  of  the  results  produced  by  RytovProp  I  have  compared 
results  produced  using  RytovProp  with  corresponding  results6  produced  using  the  split-step  wave 
optics  propagation  simulation  method,  for  the  Up-Link  Strehl  ratio.  There  were  a  total  of  66 
different  engagements  for  which  results  were  prepared.  These  engagements  differed  in  terms  of  1) 
the  satellite’s  altitude,  2)  the  zenith  angle,  3)  the  optical  wave  length,  4)  the  aperture  diameter, 
and  5)  the  optical  strength  of  turbulence.  The  satellite  altitudes  that  were  considered  were  400  km, 
1500  km,  and  that  corresponding  to  a  geo  synchronous  orbit.  The  propagation  directions  that  were 
considered  were  straight  up  (i.e.  with  a  zenith  angle  of  0  rad)  and  at  45°  to  the  zenith  direction  (i.e. 
with  a  zenith  angle  of  \  n  rad).  The  optical  wave  lengths  that  were  considered  were  A  =  0.5  /jm  and 
A  =  1.5  /jm.  The  aperture  diameters  that  were  considered  were  D  =  0.1  m  and  D  =  0.25  m.  (There 
were  wave  optics  propagation  results  developed  for  a  diameter  of  D  =  0.5  m,  but  these  results  were 
not  used  in  the  comparison  tests  being  reported  here.  The  optical  strengths  of  turbulence  that  were 
used  matched  the  HV5/7  turbulence  model — but  either  scaled  to  match  the  best  10%  integrated 
strength  observed  at  the  Starfire  Optical  Range,  or  at  its  nominal  values  (i.e.  not  scaled),  or  scaled 
to  match  the  worst  10%  integrated  strength  observed  at  the  Starfire  Optical  Range. 

For  each  engagement  three  separate  sets  of  results  were  developed,  one  set  for  a  ground  system 
having  a  tilt  tracking  servo  bandwidths  of  0  Hz  (no  tilt  tracking),  one  for  a  ground  system  with  a 
tilt  tracking  servo  bandwidth  of  3  Hz,  and  one  set  for  a  ground  system  having  a  tilt  tracking  servo 
bandwidths  of  10  Hz. 


These  results  were  prepared  by  Barry  Foucault  at  SAIC. 
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The  split-step  wave  optics  propagation  simulations  were  all  run  with  a  frame  rate  of  3,000  frames 
per  second.  These  simulations  were  run  in  sequences  of  2,800  frames,  with  the  first  700  frames  of 
each  sequence  being  discarded  as  being  “contaminated”  by  tilt  tracking  servo  start-up  transients. 
For  each  such  2,800  frame  sequence  there  was  a  single  random  number  seed  used  in  generating 
the  simulated  turbulence  patterns  that  were  used — a  random  number  seed  that  was  unique  to  that 
sequence  of  2,800  frames.  For  each  of  the  66  engagements  treated  there  were  several  such  sequences 
of  2,100  usable  frames  of  data  (i.e.  Strehl  ratio  values)  generated,  but  in  the  largest  cases  there  were 
only  20  such  sequences  generated — yielding  a  total  of  only  42,000  Strehl  ratio  values. 

It  is  inherent  in  the  way  the  split-step  wave  optics  propagation  simulation  process  is  (normally) 
used  that  there  is  in  general  very  little  change  in  the  turbulence  pattern  from  frame  to  frame.  As  a 
consequence  there  is  a  high  degree  of  correlation  between  successive  Strehl  ratio  values  produced.  In 
studying  the  frame-to-frame  correlation  of  Strehl  ratio  values  I  found  that  in  general  the  correlation 
dropped  to  a  normalized  value  of  about  one-half  with  a  separation  of  the  order  of  ten  frames — 
indicating  that  there  actually  were  no  more  than  about  4,000  degrees-of-freedom  in  any  one  of  the 
66  cases. 

The  corresponding  sets  of  RytovProp  based  results  each  had  100,000  statistically  independent  Strehl 
Ratio  values.  (It  is  inherent  in  the  RytovProp  method  that,  unless  special  steps  are  taken  to  obtain 
results  with  some  degree  of  correlation,  there  will  be  no  correlation  between  any  two  Strehl  ratio 
values  in  a  set.)  Accordingly  there  are  100,000  degrees-of-freedom  in  the  RytovProp  results  I  am 
presenting. 

Cumulative  probability  distribution  results  were  developed  for  each  of  the  66  engagement  scenarios. 
Separate  cumulative  probability  distribution  results  for  each  of  the  three  tilt  tracking  servo  band- 
widths  were  developed-  both  from  the  set  of  Strehl  ratio  values  obtained  using  the  split-step  wave 
optics  propagation  simulation  method  and  from  the  set  of  Strehl  ratio  results  obtained  using  the 
RytovProp  method.  In  general  there  appeared  to  be  good  agreement  between  the  two  sets  of  results, 
but  with  some  quite  noticeable  discrepancies  between  the  two  sets  of  results  in  the  low  probability 
range — to  some  extent  the  below  1.0%  range  but  more  so  in  the  below0.1%  range.  I  attribute  such 
discrepancies  to  the  relatively  small  number  of  degrees-of-freedom  in  the  sets  of  Strehl  ratio  values 
produced  using  the  split-step  wave  optics  propagation  simulation  method. 

To  illustrate  this  situation  in  Fig.’s  1,  2,  and  3  I  show  the  cumulative  probability  distribution  results 
for  aperture  diameters  of  D  =  0.1  m  and  D  =  0.25  m,  and  for  wave  lengths  of  A  =  0.5  /jm  and 
A  =  1.5  /jin.  The  results  shown  in  these  three  figures  were  chosen  from  the  66  different  engagement 
scenarios  studied,  were  chosen  for  presentation  here  on  the  basis  that  they  are  for  engagements 
for  which  the  values  of  the  Rytov  number — for  optical  wave  length  of  A  =  0.5  /jin  — is  about 
Vy  =  0.3  Np2,  Vy  =  0.1  Np2,  and  Vy  =  0.03  Np2.  For  Fig.  1,  i.e.  for  the  first  engagement  scenario, 
the  value  of  the  Rytov  number  for  wave  lengths  of  A  =  0.5  /jm  (and  A  =  1.5  /jin)  is  Vy  —  0.3262  Np2 
(and  Vy  —  0.0905  Np2),  while  the  effective  coherence  diameter  has  a  value  of  r0  =  0.017  m  (and 
r0  =  0.063  m),  with  the  Tyler  frequency  having  a  value  of  /T  =  56.7  Hz  (and  fT  =  18.9  Hz).  For 
Fig.  2,  i.e.  for  the  second  engagement  scenario,  the  value  of  the  Rytov  number  for  wave  lengths 
of  A  =  0.5  /jm  (and  A  =  1.5  /jm)  is  Vy  =  0.1028  Np2  (and  Vy  —  0.0285  Np2),  while  the  effective 
coherence  diameter  has  a  value  of  r0  =  0.042  m  (and  r0  =  0.158  m),  with  the  Tyler  frequency  having 
a  value  of  /T  =  28.1  Hz  (and  fT  =  9.4  Hz).  For  Fig.  3,  i.e.  for  the  third  engagement  scenario,  the 
value  of  the  Rytov  number  for  wave  lengths  of  A  =  0.5  /jm  (and  A  =  1.5  /jm)  is  Vy  =  0.0326  Np2 
(and  Vy  —  0.0091  Np2),  while  the  effective  coherence  diameter  has  a  value  of  r0  =  0.103  m  (and 
r0  =  0.0.386  m),  with  the  Tyler  frequency  having  a  value  of  /T  =  16.8  Hz  (and  /T  =  5.6  Hz). 
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X  =  0.5  |im,  D  =  0.1  m  X  =  1.5  |rm,  D  =  0.1  m 


(Fig.  la) 


(Fig.  lb) 


X  =  0.5  |im,  D  =  0.25  m 


^,=  1.5  |im,  D  =  0.25  m 


(Fig.  lc) 


(Fig.  Id) 


Figure  1.  Cumulative  Probability  Distribution  For  The  First  Engagement  Scenario 


Each  of  the  four  plots  is  for  the  combinations  of  aperture  diameter,  D,  and  optical  wave  length, 
A,  indicated  above  that  plot.  The  engagement  parameters  are  such  that  for  A  =  0.5  /Ltm  the  rel¬ 
evant  turbulence  parameters  have  values  of  Vy  =  0.3262  Np2,  r0  =  0.017  m,  and  /T  =  56.7  Hz, 
while  for  A  =  1.5  /rm  the  values  are  Vy  =  0.0905  Np2,  r0  =  0.063  m,  and  /T  =  18.9  Hz.  The 
dashed  line  curves  show  the  results  obtained  using  the  split-step  wave  optics  propagation  simula¬ 
tion  method.  The  solid  line  curves  show  the  results  obtained  using  the  RytovProp  method.  The 
red  line  curves  are  for  the  case  where  the  tilt  tracking  servo  bandwidth  was  /3dB  =  0  Hz,  while 
the  green  line  and  the  red  line  curves  are  for  /3dB  =  3  Hz  and  /3dB  =  10  Hz  respectively.  The 
horizontal  dotted  lines  indicate  the  1%  and  the  0.1%  cumulative  probability  levels. 
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X  =  0.5  |im,  D  =  0.1  m 


(Fig.  2a) 


X  =  1.5  |im,  D  =  0.1  m 


(Fig.  2b) 


X  =  0.5  |im,  D  =  0.25  m 


^,=  1.5  |im,  D  =  0.25  m 


(Fig.  2c) 


(Fig.  2d) 


Figure  2.  Cumulative  Probability  Distribution  For  The  First  Engagement  Scenario 


Each  of  the  four  plots  is  for  the  combinations  of  aperture  diameter,  D,  and  optical  wave  length, 
A,  indicated  above  that  plot.  The  engagement  parameters  are  such  that  for  A  =  0.5  /Ltm  the  rel¬ 
evant  turbulence  parameters  have  values  of  Vy  =  0.103  Np2,  r0  =  0.042  m,  and  /T  =  28.1  Hz, 
while  for  A  =  1.5  /im  the  values  are  Vy  =  0.029  Np2,  r0  =  0.158  m,  and  fT  =  9.4  Hz.  The 
dashed  line  curves  show  the  results  obtained  using  the  split-step  wave  optics  propagation  simula¬ 
tion  method.  The  solid  line  curves  show  the  results  obtained  using  the  RytovProp  method.  The 
red  line  curves  are  for  the  case  where  the  tilt  tracking  servo  bandwidth  was  /3dB  =  0  Hz,  while 
the  green  line  and  the  red  line  curves  are  for  /3dB  =  3  Hz  and  /3dB  =  10  Hz  respectively.  The 
horizontal  dotted  lines  indicate  the  1%  and  the  0.1%  cumulative  probability  levels. 
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X  =  0.5  |im,  D  =  0.1  m  X  =  1.5  |rm,  D  =  0.1  m 


(Fig.  3a) 


(Fig.  3b) 


X  =  0.5  |im,  D  =  0.25  m 


^,=  1.5  |im,  D  =  0.25  m 


(Fig.  3c) 


(Fig.  3d) 


Figure  3.  Cumulative  Probability  Distribution  For  The  First  Engagement  Scenario 


Each  of  the  four  plots  is  for  the  combinations  of  aperture  diameter,  D,  and  optical  wave  length, 
A,  indicated  above  that  plot.  The  engagement  parameters  are  such  that  for  A  =  0.5  /Ltm  the  rel¬ 
evant  turbulence  parameters  have  values  of  7 3/  =  0.033  Np2,  r0  =  0.103  m,  and  /T  =  16.8  Hz, 
while  for  A  =  1.5  /im  the  values  are  Vy  =  0.009  Np2,  r0  =  0.386  m,  and  fT  =  5.6  Hz.  The 
dashed  line  curves  show  the  results  obtained  using  the  split-step  wave  optics  propagation  simula¬ 
tion  method.  The  solid  line  curves  show  the  results  obtained  using  the  RytovProp  method.  The 
red  line  curves  are  for  the  case  where  the  tilt  tracking  servo  bandwidth  was  /3dB  =  0  Hz,  while 
the  green  line  and  the  red  line  curves  are  for  /3dB  =  3  Hz  and  /3dB  =  10  Hz  respectively.  The 
horizontal  dotted  lines  indicate  the  1%  and  the  0.1%  cumulative  probability  levels. 
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To  help  in  judging  the  significance  of  what  discrepancies  can  be  seen  in  these  figures — which  discrep¬ 
ancies  I  believe  are  almost  entirely  inherent  in  the  dashed  line  curves  and  are  due  to  the  smallness  of 
the  number  of  degrees-of-freedom  in  the  sets  of  split-step  wave  optics  propagation  simulation  Strehl 
ratio  values  used  in  forming  each  of  those  dashed  line  curves — one  needs  to  consider  not  only  the 
total  number  of  the  Strehl  ratio  values  in  each  set,  but  also  the  degree  of  correlation  of  the  successive 
Strehl  ratio  values. 

For  Fig’s,  la  and  lc  there  were  39,900  Strehl  ratio  values  used  in  forming  the  dashed  line  curves, 
while  for  Fig.’s  lb  and  Id  there  were  25,200  Strehl  ratio  values  used.  For  Fig.’s  2a  and  2c  there  were 
33,600  Strehl  ratio  values  used  in  forming  the  dashed  line  curves,  while  for  Fig.’s  2b  and  2d  there 
were  16,800  Strehl  ratio  values  used.  For  Fig.’s  3a  and  3c  there  were  23,100  Strehl  ratio  values  used 
in  forming  the  dashed  line  curves,  while  for  Fig.’s  3b  and  3d  there  were  27,300  Strehl  ratio  values 
used. 

To  provide  information  about  the  correlation  of  successive  Strehl  ratio  values  produced  by  the  split- 
step  wave  optics  propagation  method  in  Fig.’s  4,  5,  and  6  I  show  the  correlation  of  these  successive 
Strehl  ratio  values  for  the  data  sets  used  in  generating  the  curves  shown  in  Fig.’s  1,  2,  and  3 
respectively. 
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A  =  0.5  (im,  D  =  0.1  m,  Num  =  39,000 


(Fig.  4a) 


E _ I _ I  i  I  i  I  i  I  ihlihhl _ i _ I  i  I  i  I  i  1 1 1 1 1 1 1 1 1 1 1 _ i _ I  i  I  i  I  i  I  1 1 1 1 1 1 1  Ih 

io°  i  o'  io2  io3 


Frame-To-Frame  Separation 
(Fig.  4b) 


A,  =  0.5  |jm,  D  =  0.25  m,  Num  =  39,000 


(Fig.  4c) 


10°  io1  1q2  .  1t)3 

Frame-To-Frame  Separation 

(Fig.  4d) 


Figure  4.  Normalized  Correlation  Of  Successive  Strehl  Ratio  Values 


The  three  curves  in  each  plot  are  to  be  associated  with  the  similarly  colored  dashed  line  curves 
in  the  corresponding  plots  of  Fig.  1  — the  curves  based  on  the  wave  optics  propagation  method. 
The  normalized  covariance  is  the  covariance  divided  by  the  variance.  The  notation  Num  above 
each  plot  indicates  the  number  of  Strehl  ratio  values  were  used  in  producing  the  results  shown 
in  Fig.  1.  The  value  of  Num  divided  by  the  Frame-to-Frame  Separation  at  which  the  normalized 
covariance  curve  crosses  the  0.5  level  can  be  taken  as  an  estimate  of  the  number  of  degrees-of- 
freedom  in  the  data  used  forming  the  dashed  line  curves  of  Fig.  1. 
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(Fig.  5b) 


X  =  0.5  |jm,  D  =  0.25  m,  Num  =  33,600 


X  =  1 .5  |jm,  D  =  0.25  m,  Num  =  1 6,800 


(Fig.  5c) 

Figure  5.  Normalized  Correlation  Of  Successive  Strehl  Ratio  Values 


(Fig.  5d) 


The  three  curves  in  each  plot  are  to  be  associated  with  the  similarly  colored  dashed  line  curves 
in  the  corresponding  plots  of  Fig.  2  — the  curves  based  on  the  wave  optics  propagation  method. 
The  normalized  covariance  is  the  covariance  divided  by  the  variance.  The  notation  Num  above 
each  plot  indicates  the  number  of  Strehl  ratio  values  were  used  in  producing  the  results  shown 
in  Fig.  2.  The  value  of  Num  divided  by  the  Frame-to-Frame  Separation  at  which  the  normalized 
covariance  curve  crosses  the  0.5  level  can  be  taken  as  an  estimate  of  the  number  of  degrees-of- 
freedom  in  the  data  used  forming  the  dashed  line  curves  of  Fig.  2. 


-  21  - 


<D 

U 

3 

2 

& 

> 

o 

o 

<D 

N 


O 

£ 


10  10  10 

Frame-To-Frame  Separation 

(Fig.  6a) 


10 


<D 

O 

£ 

cS 

> 

O 

o 

a; 

N 


o 

£ 


10  10  10 

Frame-To-Frame  Separation 

(Fig.  6b) 
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Figure  6.  Normalized  Correlation  Of  Successive  Strehl  Ratio  Values 


The  three  curves  in  each  plot  are  to  be  associated  with  the  similarly  colored  dashed  line  curves 
in  the  corresponding  plots  of  Fig.  3  — the  curves  based  on  the  wave  optics  propagation  method. 
The  normalized  covariance  is  the  covariance  divided  by  the  variance.  The  notation  Num  above 
each  plot  indicates  the  number  of  Strehl  ratio  values  were  used  in  producing  the  results  shown 
in  Fig.  3.  The  value  of  Num  divided  by  the  Frame-to-Frame  Separation  at  which  the  normalized 
covariance  curve  crosses  the  0.5  level  can  be  taken  as  an  estimate  of  the  number  of  degrees-of- 
freedom  in  the  data  used  forming  the  dashed  line  curves  of  Fig.  3. 
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Over  all  the  agreement  between  the  RytovProp  results  and  the  split-step  wave  optics  propagation 
results  shown  in  Fig.’s  1,  2,  and  3  is  quite  good,  especially  when  allowance  is  made  for  the  smallness 
of  the  number  of  degrees-of-freedom  of  the  data  sets  used  to  produce  the  wave  optics  propagation 
probability  distribution  curves  (the  dashed  line  curves).  I  take  this  as  tending  to  validate  the 
RytovProp  method. 
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'  A  covariance  matrix  can  not  have  negative  eigen-values 


